library(tidyverse)
library(rio)
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(geobr)
library(sf)

bd <- import("G:/.shortcut-targets-by-id/1WF0Sd5YrrhwmtJxmTJQyGG_fVdGtSIWt/Editoria Dados/2025 n 2   n 3 NOVA REMESSA/2025 n 3/Artigo X/Base de dados - Federalismo e desigualdades horizontais em educação.xlsx")

# Tabela 1 ---------------------------------------------------------------------

names(bd)

summary(bd)

bd |> 
  group_by(cwz5) |>
  summarise(n = n()) |> 
  na.omit(czw5) |> 
  mutate(total = sum(n),
         percent = n/total * 100)

bd |> 
  group_by(cwz5) |>
  summarise(n = n()) |> 
  na.omit() |> 
  mutate(total = sum(n),
         percent = n/total * 100,
         cumulative_percent = cumsum(percent))

# Tabela 2 ---------------------------------------------------------------------
# Incluir região 

df <- read_csv("br_bd_diretorios_brasil_municipio.csv.gz")

regiao <- df %>%
  select(id_municipio, sigla_uf, nome_regiao, nome)


bd <- left_join(bd, regiao, by = c("mun" = "id_municipio"))


names(bd)



# Calcular o número total de municípios em cada cluster
total_por_cluster <- bd %>%
  filter(!is.na(cwz5)) %>%  # Filtrar NAs no cluster
  group_by(cwz5) %>%
  summarise(total_cluster = n())

# Calcular o número de municípios em cada cluster por região
tabela <- bd %>%
  filter(!is.na(cwz5), !is.na(nome_regiao)) %>%  # Filtrar NAs no cluster e região
  group_by(nome_regiao, cwz5) %>%
  summarise(n = n(), .groups = 'drop') %>%
  left_join(total_por_cluster, by = "cwz5") %>%
  mutate(percent = n / total_cluster * 100) %>%
  select(-n, -total_cluster) %>%
  pivot_wider(names_from = cwz5, values_from = percent, values_fill = 0) %>%
  mutate(Total = rowSums(across(where(is.numeric)), na.rm = TRUE))





# Calcular o total geral para cada cluster
total_geral <- tabela %>%
  summarise(across(-nome_regiao, ~ sum(.x, na.rm = TRUE))) %>%
  mutate(nome_regiao = "Todos", Total = 100)

# Adicionar a linha de totais gerais à tabela
tabela <- bind_rows(tabela, total_geral)

tabela


# Tabela 3 ---------------------------------------------------------------------

# Calcular o número total de municípios em cada região


total_por_regiao <- bd %>%
  filter(!is.na(nome_regiao)) %>%  # Filtrar NAs na região
  group_by(nome_regiao) %>%
  summarise(total_regiao = n())

# Calcular o número de municípios em cada região por cluster
tabela_invertida <- bd %>%
  filter(!is.na(cwz5), !is.na(nome_regiao)) %>%  # Filtrar NAs no cluster e região
  group_by(nome_regiao, cwz5) %>%
  summarise(n = n(), .groups = 'drop') %>%
  left_join(total_por_regiao, by = "nome_regiao") %>%
  mutate(percent = n / total_regiao * 100) %>%
  select(-n, -total_regiao) %>%
  pivot_wider(names_from = cwz5, values_from = percent, values_fill = 0) %>%
  mutate(Total = rowSums(across(where(is.numeric)), na.rm = TRUE))

# Calcular o total geral para cada região
total_geral_invertido <- tabela_invertida %>%
  summarise(across(-nome_regiao, ~ sum(.x, na.rm = TRUE))) %>%
  mutate(nome_regiao = "Todos", Total = 100)

# Adicionar a linha de totais gerais à tabela
tabela_invertida <- bind_rows(tabela_invertida, total_geral_invertido)

tabela_invertida


### Teste adicional 

# Calcular o total de municípios por região que estão em algum cluster
total_por_regiao <- bd %>%
  filter(!is.na(nome_regiao), !is.na(cwz5)) %>%  # Certifique-se de considerar apenas municípios com cwz5 não NA
  group_by(nome_regiao) %>%
  summarise(total_regiao = n(), .groups = 'drop')

# Calcular o número de municípios em cada região por cluster, ajustando a base de dados inicial
tabela_invertida <- bd %>%
  filter(!is.na(cwz5), !is.na(nome_regiao)) %>%
  group_by(nome_regiao, cwz5) %>%
  summarise(n = n(), .groups = 'drop') %>%
  left_join(total_por_regiao, by = "nome_regiao") %>%
  mutate(percent = n / total_regiao * 100) %>%
  pivot_wider(names_from = cwz5, values_from = percent, values_fill = list(percent = 0)) %>%
  mutate(Total = rowSums(across(where(is.numeric)), na.rm = TRUE))

# Visualizar a tabela para confirmar que os totais agora estão corretos
print(tabela_invertida)


###

# Calcular o total de municípios por região que estão em algum cluster
total_por_regiao <- bd %>%
  filter(!is.na(nome_regiao), !is.na(cwz5)) %>%
  group_by(nome_regiao) %>%
  summarise(total_regiao = n(), .groups = 'drop')

# Calcular o número de municípios em cada região por cluster e formatar a tabela
tabela_invertida <- bd %>%
  filter(!is.na(cwz5), !is.na(nome_regiao)) %>%
  group_by(nome_regiao, cwz5) %>%
  summarise(n = n(), .groups = 'drop') %>%
  left_join(total_por_regiao, by = "nome_regiao") %>%
  mutate(percent = n / total_regiao * 100) %>%
  select(nome_regiao, cwz5, percent) %>%
  pivot_wider(names_from = cwz5, values_from = percent, values_fill = list(percent = 0), names_prefix = "") %>%
  mutate(Total = rowSums(across(where(is.numeric)), na.rm = TRUE))

# Visualizar a tabela final
print(tabela_invertida)


### Figura 1 -------------------------------------------------------------------



# download municipality geometries
muni_sf <- geobr::read_municipality(year = 2010,
                                    showProgress = FALSE)
#> Using year/date 2010

# merge data
muni_sf$code_muni <- as.double(muni_sf$code_muni)
esg_sf <- left_join(muni_sf, bd, by = c('code_muni' = "mun"))

# Definir a paleta de cores
colors <- c("1" = "yellow", 
            "2" = "orange", 
            "3" = "darkorange", 
            "4" = "lightcoral", 
            "5" = "darkred")

# Criar o mapa
ggplot() +
  geom_sf(data = esg_sf, aes(fill = as.factor(cwz5)), color = NA) +
  scale_fill_manual(values = colors, na.value = "white") +
  labs(fill = "Cluster") +
  theme_minimal() +
  theme(legend.position = "right")



# Tabela 4 ---------------------------------------------------------------------
# verificando apenas algumas das linhas

bd |> 
  group_by(cwz5) |>
  summarise(media_matriculas = mean(NU_MATRICULAS, na.rm = T),
            media_prova_brasil_ai = mean(prova_brasil_ai, na.rm = T),
            rotatividade = mean(cl_rotati_alta, na.rm = T)) |> 
  na.omit(czw5)



# Tabela 5 ---------------------------------------------------------------------

bd |> 
  group_by(cwz5) |>
  summarise(media_efic_ai = mean(vrste_ai, na.rm = T),
            media_efic_af = mean(vrste_af, na.rm = T)) |> 
  na.omit(czw5)


# Tabela 6 ---------------------------------------------------------------------

bd |> 
  group_by(efi_ai) |>
  summarise(n = n())

bd |> 
  group_by(efi_af) |>
  summarise(n = n())


# Tabela 7 ---------------------------------------------------------------------
# verificando apenas algumas das linhas

# gasto médio por aluno

table(bd$efi_ai)

bd |> 
  group_by(efi_ai) |>
  summarise(media_gasto_aluno = mean(gmpa2014, na.rm = T))

bd <- bd %>%
  mutate(efi_ai_grupo = case_when(
    efi_ai %in% c("Eficientes", "Ineficiencia fraca") ~ "Eficientes/Ineficiencia fraca",
    efi_ai == "Ineficiencia forte" ~ "Ineficiencia forte",
    efi_ai == "Ineficiencia moderada" ~ "Ineficiencia moderada",
    TRUE ~ NA_character_
  ))


# Calcular a média do gasto por aluno para os grupos desejados
media_gasto <- bd %>%
  filter(efi_ai_grupo %in% c("Eficientes/Ineficiencia fraca", "Ineficiencia forte")) %>%
  group_by(efi_ai_grupo) %>%
  summarise(media_gasto_aluno = mean(gmpa2014, na.rm = TRUE))

# Ver o resultado
print(media_gasto)

# nível socioeconômico
bd |> 
  group_by(efi_ai) |>
  summarise(nivel_socioeconomico = mean(INSEVALORABSOLUTO, na.rm = T))

# Calcular a média do nível socioeconômico para os grupos desejados
media_socio <- bd %>%
  filter(efi_ai_grupo %in% c("Eficientes/Ineficiencia fraca", "Ineficiencia forte")) %>%
  group_by(efi_ai_grupo) %>%
  summarise(media_nivel_socio = mean(INSEVALORABSOLUTO, na.rm = TRUE))

# Ver o resultado
print(media_socio)



# infraestrutura
bd |> 
  group_by(efi_ai) |>
  summarise(infraestrutura = mean(cl_infra, na.rm = T))

# Calcular a média da infraestrutura para os grupos desejados
media_infra <- bd %>%
  filter(efi_ai_grupo %in% c("Eficientes/Ineficiencia fraca", "Ineficiencia forte")) %>%
  group_by(efi_ai_grupo) %>%
  summarise(media_infraestrutura = mean(cl_infra, na.rm = TRUE))

# Ver o resultado
print(media_infra)


# Tabela 8 ---------------------------------------------------------------------

bd |> 
  group_by(efi_ai) |>
  summarise(n = n(),
            media_prova_br = mean(prova_brasil_ai, na.rm = T))


bd |> 
  group_by(efi_af) |>
  summarise(n = n(),
            nivel_socioeconomico = mean(INSEVALORABSOLUTO, na.rm = T))


# Tabela 9 ---------------------------------------------------------------------

library(texreg)

# Base 1: Definir o cluster 3 como referência
bd_3_ref <- bd |> 
  mutate(cwz5 = relevel(as.factor(cwz5), ref = "3"))

# Ajustar o modelo para vrste_ai com cluster 3 como referência
modelo_ai <- lm(vrste_ai ~ cwz5, bd_3_ref)
summary(modelo_ai)

# Base 2: Definir o cluster 4 como referência
bd_4_ref <- bd |> 
  mutate(cwz5 = relevel(as.factor(cwz5), ref = "4"))

# Ajustar o modelo para vrste_af com cluster 4 como referência
modelo_af <- lm(vrste_af ~ cwz5, bd_4_ref)
summary(modelo_af)

# Exibir os resultados dos modelos em conjunto
screenreg(list(modelo_ai, modelo_af))


# RASCUNHO ---------------------------------------------------------------------

bd |> 
  na.omit(cwz5) |> 
  count(nome_regiao)

# Tabela 3 (Carlos)
tabela_3 <- bd |> 
  na.omit(cwz5) |> 
  group_by(nome_regiao, cwz5) |> 
  summarise(n_municipios = n()) |> 
  mutate(total_municipios = case_when(nome_regiao == "Centro-Oeste" ~ 158,
                                      nome_regiao == "Nordeste" ~ 1185,
                                      nome_regiao == "Norte" ~ 155,
                                      nome_regiao == "Sudeste" ~ 642,
                                      nome_regiao == "Sul" ~ 297),
         percent = n_municipios / total_municipios * 100)
